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Abstract 

I In this paper we construct a third order method for solving additively split autonomous stiff 

■ systems of ordinary differential equations. The constructed additive method is L-stable with respect 

(-H I to the implicit part and allows to use an arbitrary approximation of the Jacobian matrix. Automatic 

stepsize selection based on local error and stability control are performed. The estimations for 
error and stability control have been obtained without significant additional computational costs. 
Numerical experiments show reliability and efficiency of the implemented integration algorithm. 



> ' 1 Introduction 

Q\ . Spatial discretization of continuum mechanics problems in partial differential equations by finite 

CN I difference or finite element methods results in the Cauchy problem for the system of ordinary differential 

CN ' equations with an additively split right hand side function of the form: 
O " 

o 



y' = ^{i,y) + 9ii,y), y{to)=yo, to<t<tk, 



where (p{t, y) is a non-symmetrical term obtained from discretization of the first-order differential 
^ ' operator, g{t, y) is a symmetrical term obtained from discretization of the second-order differential 

operator, i is a independent variable. It is assumed that in the problem the vector-function g is a stiff 
term and (f is a non-stiff term. 

Explicit Runge-Kutta methods have a bounded stability region and are suitable for non-stiff and 
mildly stiff problems only. L-stable methods are usually used for solving stiff problems. In the case of 
large-scale problems overall computational costs of L-stable methods are almost completely dominated 
by evaluations and inversions of the Jacobian matrix of a right hand side vector function. Overall 
computational costs can be significantly reduced by re-using the same Jacobian matrix over several 
integration steps (freezing the Jacobian). 

Freezing the Jacobian in iterative methods has effect on convergence speed of an iterative process 
only and doesn't lead to loss of accuracy. So, this approach is extensively used for implementation 
of these methods. For Rosenbrock type methods and their modifications [4] an approximation of the 
Jacobian matrix can lead to decreasing a consistency order. 

The system y' = f{t,y) can be written in the form y' = [f{t,y) — By] + By, where B is 
some approximation of the Jacobian matrix. Assume that stiffness is fully concentrated in the term 
g{t,y) = By, then the expression ip{t,y) = f{t,y) — By can be interpreted as the non-stiff term [2,7]. 
If the Cauchy problem is considered in the form y' = [f{t, y) — By] + By under construction of additive 
methods, then an arbitrary approximation of the Jacobian matrix can be used without decreasing the 
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order of these methods. Additive methods constructed in this way allow both analytical and numerical 
computations of the Jacobian matrix. Note that the approximation of the Jacobian by a diagonal 
matrix is suitable for some mildly stiff problems. 

In this paper we construct a six-stage third order additive method that allows to use different 
kinds of approximation of the Jacobian matrix. The estimations of the error and the maximum absolute 
eigenvalue of the Jacobian matrix have been obtained without significant additional computational 
costs. Indeed, the error estimation has been obtained on the base of an embedded additive method 
and the maximum absolute eigenvalue estimation has been obtained by a power method using only 
two additional computations of ^{y). These estimations are used for error and stability control 
correspondingly. Numerical experiments are performed showing the reliability and efficiency of the 
constructed method. 



2 A numerical scheme for autonomous problems 

Consider the Cauchy problem for an autonomous system of ordinary differential equations 

y' = ^{y) + 9{y), y{io) =yo, to <t < tk, 



where y, if and g are A^- dimensional smooth vector-functions, t is an independent variable. In the 
following, we assume that g is a stiff term and is a non-stiff term. Consider a six-stage numerical 
scheme for solving (1): 



2/n+l =yn + ^Piki, 
i=l 

h = hip{yn), 

Dnk2 = h[ip{yn) +g{yn)], 

Dnh = h, 

3 3 

Dnki = h(p{yn + ^ f^ijkj) + hg{yn + ^ aijkj), (2) 

i=i 3=1 
Dnk5 = k4 + 'jks, 

5 

ke = hip{yn + ^ (3&jkj), 



where -D„ = E — ahg'^, g'^ = dg{yn)/dy is the Jacobian matrix of the function g{y), E is the identity 
matrix, fcj, 1 < i < 6, are stages, a,pi,a4j, l34^j, Pqj,^ are coefficients that have effect on accuracy and 
stability properties of the scheme (2). 



3 The third order conditions 



The Taylor series expansion of the approximate solution up to terms in has the form 



Vn+l =yn+ {Pl + P2 + P3 + P4 + (7 + 1)P5 + Pejhip + (p2 + PS + P4 + (T + ^)P5)hg+ 
+ ((/341 + /342 + /343)(P4 + Ps) + {Pel + P%2 + /^eS + ^64 + (t + 1)^65)P6) 

+ ((/342 + /343)(P4 + P5) + (/562 + /^GS + /364 + (t + l)/365)p6) h^v'g+ 

+ [a{P2 + 2p3 +PA + (37 + 2)p5) + (041 + 042 + a43)(P4 + P5)] h^g''f+ 
+ [a(p2 + 2p3 + P4 + (37 + 2)^5) + ("42 + a43)(p4 + P5)] q' 9+ 

+ 0.5[(/34i + I3a2 + PAZ?{PA + P^) + (/561 + /962 + /^es + /364 + (7 + l)/365)^P6]/lVV^+ 
+ 0.5 [(/342 + /?43)'(P4 +P5) + (/362 + /^es + /?64 + (7 + I) Pfi^f PQ\h\" + 

+ [(/342 + /343)(/941 + (3a2 + Pa3)(.PA + P5) + (/?61 + ^62 + ^63 + p6A + 
+ (7 + l)/365) (/362 + /363 + PgA + (7 + 1)/^65)P6] /^VVs + iPil + p42 + 

+ /343)(/964 + /365)P6/iW + (/342 + pA3){Pe4 + Pe5)peh^v'^ 9+ 



+ 



+ 



((/342 + 2/343)(P4 + Ps) + {P62 + 2/303 + ^64 + (87 + 2)/365)p6) " 
+ (a4i + a42 + a43)(/564 + /365)P6 h^v'g'f + a(^(/342 + 2/343)(P4 + P5) + 

+ {p62 + 2/363 + /364 + (37 + 2)/365)P6) + ("42 + 043) (/364 + p65)P6 h^v'g'g+ 

+ 0.5(q4i + a42 + a43)^(P4 + P5)h^9"v>'^ + 0.5(a42 + a43)^(P4 + P5)h^9"g^+ 

+ {aA2 + 043) ("41 + "42 + 043) (P4 + P5)h^g"^g + a(/341 + Pa2 + /343)(P4+ 

+ 2p5)h^g'ip'ip + a(/342 + /343)(P4 + 2p5)h^g'ip'g + a[a{p2 + 3p3 + pa+ 

+ (67 + 3)^5) + (041 + 2^42 + 3q;43)?'4 + (2a4i + 2^42 + 3a43)P5] h^g''^'P+ 

+ a[a{p2 + 3p3 +P4 + (67 + 3)p5) + (2042 + 3043)^4 + (2a42 + 3a43);j5] h^g'^g+ 

+ Oih% 



where the corresponding elementary differentials are evaluated at i/n- 

The Taylor series expansion of the exact solution up to third order terms is 



/j2 ^3 
y(Wi) = yitn) + h{<p + g) + —{ip'ip + if'g + g'ip + g'g) + y (¥''V'+ 



+ ip"g'^ + 2ip"ipg + + ip'^g + ip'g'ip + ip'g'g + g"ip^ + g"g^+ 

+ '^g"'pg + ffW + g'f'g + g'^'P + g'^g) + o{h'^), 



i2 



l',J2 



J> „2 



(3) 



where the corresponding elementary differentials are evaluated at y{tn). 

Comparing the successive terms in the Taylor series expansion of the approximate and the exact 
solutions up to third order terms under the assumption ?/„ = y{tn) we have the system of nonlinear 
algebraic equations. Its solving results in the relation Pai = olai = Pei = and the third order 



conditions of the scheme (2) take form: 



1) P2 + P3 + P4 + (7 + 1)P5 = 1, 

2) (/?41 + P42 + /343) (P4 + Ps) + (/?61 + P62 + p63 + /?64 + 
+ (7 + 1)/365)P6 = 0.5, 

3) a{p2 + 2p3 +P4 + (37 + 2)^5) + (a4i + 042 + a43)(p4 +P5) = 0.5, 

4) (/341 + ^42 + PisfiPA + P5) + (/?61 + /?62 + PgS + /364+ 
+ (7+l)/365)'p6 = l/3, 

5) (/?41 + /342 + /343)(/364 + /365)P6 = 1/6, 

6) a [{(342 + 2/343) (P4 + P5) + (/362 + 2/363 + /364 + (37 + 2)/365)p6] + (4) 

+ (041 + 042 + 043) (/364 + /365)P6 = 1/6, 

7) (a4i + a42 + a43)^(p4 +P5) = 1/3, 

8) a(/34i + /342 + Pa3){Pa + 2p^) = 1/6, 

9) a[a(p2 + 3p3+?'4 + (67+3)p5) + (a4i + 2a42+3a43)P4+ 

+ (2041 + 3042+4043)^5] = 1/6, 
10) 041 = I3ii = /361 = 0, pi = -pQ. 

4 Stability analysis 

The hnear stabihty analysis of the additive scheme (2) is based on the scalar model equation 

y' = Xiy + Aay, y(0) =2/0, ^ > 0, Re(Ai) < 0, Re(A2) < 0, | Re(Ai)| « | Re(A2)|, (5) 

where the free parameters Ai, A2 can be interpreted as some eigenvalues of the Jacobian matrices of 
the functions ip (the non-stiff term) and g (the stiff term) correspondingly. 

Application of the scheme (2) for numerical solving the equation (5) yields 

Vn+i = R{x,z)yn, 

where x = Ai/i, z = A2/1 and R{x,z) is a stability function (its analytical expression is omitted here 
for brevity). 

The necessary condition of L-stability of the additive scheme (2) with respect to the stiff term 
has the form: 



lim R{x, z) = 0. 
2— >— 00 

It is satisfied if the following two conditions hold: 

a^iPi +P&) + ((042 - a)^64 - af3Q2)pG = 0, 
a{a - P2) + (042 - a)p4 = 0. 



(6) 



Solving the system (4), (6). In the following, we assume that Y^j=i '^ij = 1 1 0!42 = Q^, Pa2 = a. 

3 

The first relation ensures that g{yn + Yl ^^Aj^j) approximate g{y{tn+i)) in the fourth stage and the 

i=i _ 

other ones improve stability properties of the intermediate numerical formulas. 



Let us denote 



Pi := /364 + P65, h ■■= Pes + P64 + (7 + l)/365, 

Ps ■■= a(2/363 + /?64 + (37 + 2)/?65) + ^64 + /?65, /34 := a + /343- 

Then after obvious simplifications the system (4), (6) takes the form 

1) P2 +P3 + 7P5 = 2/3, 

2) api{pi + 2p5) = 1/6, 

3) /34/3iP6 = 1/6, 

4) Pa/3 + p2P(i = 0.5, 

5) /3|/3 + /3|p6 = 1/3, (7) 

6) a(2/34 - a)/3 + /^sPe = 1/6, 

7) P4 +P5 = 1/3, 

8) a{p2 + 2p3 + P4 + (37 + 2)P5) = 1/6, 

9) a(a(p2 + 3p3 +P4 + (67 + 3)^5) + (3 - a)p4 + (4 - 0)^5) = 1/6, 

10) a4i = P41 = Pel = P62 = 0, 042 = /?42 = a, 043 = 1 - a, p2 = a, pi = -pQ. 

Let us consider the equations (1) and (7) - (10). Multiply (1) by 3a and subtract the result from (8). 
Then divide (9) by a and subtract (1) multiplied by 6a from it. As the result we obtain: 

P2 +P3 +7P5 = 2/3, 
PA+Vb = 1/3, 

a(-2p2-P3+P4 + 2p5) = l/6-2a, (8) 
a(-5p2 - 3p3 +PA + 3p5) + (3 - a)pA + (4 - a)p^ = (6a)"^ - 4a, 
P2 = a. 

Prom the second equation of (8) we have p^ = 1/3 — p^. Substituting this relation and the fifth relation 
of (8) to the first tree equations we have 

P3 + 7P5 = -a + 2/3, 

a{-P3 + P5) = 2a2 - 7a/3 + 1/6, (9) 
a{-3p3 + 2p5) +P5 = Sa^ - 4a - 1 + (6a)~\ 

Prom the second equation of (9) we have 



P3 =P5 



(12a2 - 14a + l)(6a)~^ (10) 



It follows from (10) and the third equation of (9) that = (6a^ - 18a^ + 9a — l)/(6a^ — 6a). 
Substituting p^ to (10) we obtain p^ = (3a^ — 4a + 3)/(3 — 3a). Substituting p^ to (10), we have 
P3 = (3a^ — 4a + 3)/(3 — 3a). Then we substitute p^ and ps to the first equation of (9). As the result we 
have 7 = 2a(a + l)/(6a^ — 18a^ + 9a — 1). It follows from the second equation of (8) that pa = (6a^ — 
20a2 + lla-l)/(6a-6a2). From the second equation of (7) we obtain Pa = (a-l)/(6a^-16a2 + 7a-l). 
The fourth and fifth equations of (7) can be written in the form P2P6 = {3 — 2Pa)/6, P2P6 = {1 — Pa)/'^- 
Dividing the second equation by the first one results in P2 = {2 - 2/5|)/(3 - 2Pa). It follows from here 



and the fourth equation of (7) that pq = (0.5 — l3i/2>)/(32- Prom the third and sixth equations of (7) 
we obtain j3i = {QISa^Pq)'^ and (3^ = [1/6 — a(2/?4 — /342)/3]/p6- 

Then we express the coefficients /Sga, /364 and /Sgs of the scheme (2) in terms of the auxihary 
parameters /3i, /?2, /Ss and /34, that is 

/564 + /365 =/3l, /363 + 7/365 =/32-/3i, a(/?63-/?65) =a(3/?2-2/?i)+/3i-/33. (11) 

Multiplying the second equation of (11) by a and subtracting it from the third one we obtain /Jgs = [a(/9i- 
2/32) + /Ss — (07 + a)- It follows from here and the second equation of (11) that /^ea =1^2 — ^1 — ll^&bi 
and from the first relation we obtain /364 = (3i — 065- 

As the result all the coefficients of the scheme (2) are expressed in terms of one free parameter a. 
The coefficient a can be found from following considerations. Let if{y) = 0, that is consider the system 
of the form y' = g{y) instead of (1). In this case the local error 5n+i at point tn+i can be represented 
in the form 

5n+i = h\c^g'^g + C2g"g'g'' + c^g'g"g'' + c^g"'g^) + 0{h'), 

where Cj, 1 < i < 4, are expressed in terms of the coefficients of the scheme (2) (their expressions are 
omitted here for brevity). 

The system y' = g{y) is stiff, that is the function g{y) satisfies the Lipschitz condition with a large 
constant. Therefore the term cih^g'^g makes the largest contribution to the local error. Choose ci = 
for minimizing the local error, then we have 

24a^ - 96a^ + 72a^ - 16a + 1 = 0. (12) 
Now, the coefficients of the L -stable third order scheme (2) can be computed by the following formulas: 

a4i = f3ii = /^ei = /362 = 0, 042 = /342 = a, 043 = 1 - a, p2 = a, 

7 = 2a{a + l)/(6a^ - 18a^ + 9a - 1), ps = (a^ - 4a/3 + 1)/(1 - a), 
P4 = (6a3 - 20a2 + 11a - l)/(6a - Ga"^), 
P5 = (6a^ - ISa^ + 9a - l)/(6a2 - 6a), 

/J4 = (a - l)/(6a^ - 16a2 + 7a - 1), I3^^ = p^-a, (13) 

/32 = (1 - /?|)/(l-5 - /34), P6 = (0.5 - /34/3)//32,Pi = -Pe, 

Pi = {QPiP6r\ 133 = [1/6 - a(2/34 - /342)/3]/p6, 
/365 = HPi - 2P2) + P3- /3i]/(a7 + a), 

Pes = P2- Pi - 7/^65, Pm = Pi - P&b, 



where the coefficient a is determined from the equation (12). 

This equation has the four real roots ai = 0.10643879214266, a2 =0.22042841025921, 
as =0.57281606248213 and 04 =3.1003167351160. The numerical experiments that we have done show 
that the root is the most suitable. Therefore computational results will be given for a = 03. The 



corresponding coefficients of the L -stable third order scheme (2) take the form 



a 


= +0.57281606248213, 


Pi = 




0.48695861160293, 


P2 


= +0.57281606248213, 


P3 


= +1.32112526220103, 


P4 = 




0.09105090402502, 


P5 


= +0.42438423735836, 


P6 


= +0.48695861160293, 


041 = 


0: 


1 


"42 


= +0.57281606248213, 


"43 


= +0.42718393751787, 


PAI = 


0: 


) 




= +0.57281606248213, 


/343 


= -0.18882050162852, 


Pel = 


0: 


1 


fim 


= 0, 


/?63 


= +2.51499368618962, 






0.022405291307077, 






/365 


= +0.91371881359685, 


7 = 




2.891895009239397. 







(14) 



5 Local error estimation 

For the error estimation we construct the embedded method of second order of the form: 

4 ^ 
Vn+l, 2 = ?/n + ^ riki + rs/cs , 
1=1 

ki = hip{yn), 

Dnk2 = h{ip{yn) +g{yn)), 

Dnk-s = k2, (15) 
3 3 

Dnki = hip{yn + ^ Pijkj) + hg{yn + ^ oiAjkj), 
j=i i=i 

Dnk^ = ^4, 

where the coefficients rj, 1 < i < 5, should be determined, and parameters a, a4j,/34j are given by (13) 
or (14). Note that there is not a sixth stage in (15) and there is not 7/03 in the fifth stage as opposed 
to (2). 

The Taylor series expansion of the approximate solution computed by the scheme (15) up to terms 
in has the form 

Vn+i, 2 = yn + (r-1 + r2 + rs + r4 + r5)/iv? + (r2 + rs + r4 + rr^)hg+ 

+ (a(r2 + 2r3 + r4 + 2r^) + r4 + r^)h^g'ip + (a(r2 + 2r3 + r4 + 2r^) + 
+ r4 + r^)h^g'g + j3i{ri + r^)h?ip'^ + p^ir^ + r^)h^^'g + 0{h^), 

where the elementary differentials are evaluated at y^. Comparing successive terms in the Taylor series 
expansion of the approximate and the exact solutions up to second order terms under the assumption 
Un = y{tn) we obtain the second order conditions of the scheme (15): 

1) ri + r2 + r3 + r4 + rs = 1, 

2) r2 + r3 + r4 + rs = 1, (16) 

3) /34(r4 + r5) =0.5, 

4) a(r2 + 2r3 + r4 + 2r5) + r4 + rs = 0.5, 

where /34 is determined by (13). Note that it follows from the ffist two equations of (16) that ri = 0. 
Now we analyze the stability of the scheme (15). Its application for numerical solving the equation (5) 
yields 

yn+1,2 = R2{x,z) yn,2 , 



where x = Xih, z = A2/1 and the stability function R2{x, z) has the form 



R2{x, z) = [a^{a — r2)z^ — a?{r2 — r^)xz'^ — a(4a^ — a(3r2 + + 2r4) + r4)z'^+ 
+ a^r^^x^z^ + a(a(3r2 + ra + r4 — rs) - r4(/34 + 

+ (6a^ - a(3r2 + 2r3 + 3r4 + 2rr^) + r4 + r^)z^ - a[a{r4 + rs) + r4/34)x^z+ 
+ (-a(3r2 + 2r3 + 3r4 + 2r5) + {r^ + r5)(/?4 + l))xz+ 

+ (-4a + r2 + rs + r4 + rs)^ + P^ir^ + r5)x^ + (r2 + rs + r4 + r5)a; + l]/(l-a2;)^. 

The necessary condition of L-stabihty of the additive scheme (15) with respect to the stiff term has 
the form: 

hm R2{x, z) = 0. 
2— >— 00 

It is satisfied if r2 = a. 

Now we consider the system (16). Dividing the third equation by (j^ and subtracting the result from 
second one we obtain = 1 — a — O.S/?^^^. Expressing in terms of from the third equation (16) 
and substituting it to the fourth equation of (16) we have r4 = 0.5(1 - /?4)(a/34)"^ + 2 - a, rg = 0.5(a- 
1 + p4){a(5i)-^ - 2 + a. 

As the result we have all the coefficients of the L stable embedded method (15) of second order. 
For the coefficients (14) we obtain 

n = 0, r2 = +0.57281606248213, rg = -0.87491444843356, 

r4 = +2.82745609901376, = -1.52535771306233. 

The embedded method (15) requires, at each integration step, only one additional backward 
substitution steps of Gauss elimination method and doesn't require additional computations of right 
hand side, evaluations and inversions of the Jacobian matrix. In the case of large-scale problems overall 
computational costs of the method (15) are almost completely dominated by evaluations and inversions 
of the Jacobian matrix. So, we obtain the error estimation based on the embedded method (15) without 
significant additional computational costs. 

Let us denote the error estimation by 

Wn ~ J/n, 2I 

ervn = max — — ; — — — r-, 

i<i<N Atoli + Rtoli\yl\ 

where Atoli and RtoU are the desired tolerances prescribed by the user. If err„ < 1, then the computed 
step is accepted, else the step is rejected and computations are repeated. When RtoU = 0, the absolute 
error is controlled on the z-th component of the solution with the desired tolerance Atoli. If Atoli = 
then the relative error is controlled on the i-th component with the tolerance Rtoli. 



6 Stability control and stepsize selection 



In the additive method (2) for solving (1) the non-stiff term tp is treated by the tree-stage exphcit Runge 
Kutta method (the exphcit part), and the stiff term g is treated by the L-stable (4, 2)-method [10-12] 
(the imphcit part). In the general case there is no guarantee that the function 99(2/) = f{y) — By is the 
non-stiff term in reducing y' = f{y) to y' = [f{y) — By] + By. If some stiffness is in ip{y) = f{y) — By 
(i.e. stiffness leakage phenomenon occurs) then the additional stability control of the explicit part of 
the scheme (2) can increase efficiency of computations for many problems. In some cases it has no a 
significant effect on the efficiency of the integration algorithm because of the good stability properties 
of the scheme (2). Therefore the choice of using or not using the additional stability control of the 
explicit part is given to the end-user. 

We perform the stability control of the explicit part of the scheme (2) by analogy with [9]. For 
additive methods in opposite to explicit Runge Kutta methods it isn't possible to use previously 
computed stages because of peculiarity of the problem (1). Therefore instead of using the stages 
h, 1 < i < 6, of (2) we consider the additional stages di, ^2 of the form: 

di = hipiyn + a2iki), c?2 = hip{yn + asih + a32di). 

Denote ^p{y) = Ay + b, where A and b are matrix and vector with constant coefficients correspondingly, 
then we have 

ki = h{Ayn + b), di = ki + a2ihAki , d2 = ki + {aai + a32)hAki + a2ia32h^A^ki. 

Assuming 021 = 031 + 032 we obtain 

d2 — di = a2ia32h^A'^ki, di — ki = a2ihAki. 

The maximum absolute eigenvalue Vn = h\Xn max\ of the matrix hA can be approximated using the 
power method by the following formula: 

-1, 14-41 



1032 I max 



i<i<N \d\-k\\ 

then the stability control can be made by f„ < 2, where number 2 is an approximate length of the 
stability interval of the tree-stage explicit Runge Kutta method. 

In the general case this estimation is quite crude because of small number of iterations of the 
power method and the nonlinearity of the function ip{y). Therefore the stability control is used for 
limiting the stepsize growing only. 

Let the approximate solution y„ is computed with the stepsize hn- For the stepsize selection we 
use errn = 0{h^). The stepsize hacc predicted by accuracy we compute by the formula: hacc = Qihn, 
where qi is a root of the equation qferVn = 1- In view of Vn = 0{hn), the stepsize hst predicted by 
stability is computed by h^t = 92^n, where 52 is a root of the equation q2Vn = 2. Then the stepsize 
predicted by accuracy and stability is selected by the formula: 

hn+i = max[hn,inm{hacc,hst)]- 



The stability control of the explicit part of the scheme (2) requires, at each integration step, two 
additional computations of (p{y). These computational costs are negligible for large-scale problems, but 
if you are sure that all stiffness is in g{y) then you can take off stability control to save computational 
costs. 



7 Numerical experiments 



Further, the numerical code based on the additive method (2) (with error and stabihty control as 
well as with diagonal Jacobian approximation) is called AS0DE3 (the Additive Solver of Ordinary 
Differential Equations) . 

The test problems given below have been reduced to the form y' = {f{y) — By) + By. All numerical 
computations have been performed in double precision arithmetic on IBM PC Athlon(tm) XP 2000+ 
with the desired tolerances of the error Atol = Rtol = Tol = 10~™, m = 2,4. The scheme (2) is of 
third order, therefore it is unreasonable to do numerical computations with higher tolerance. 

The following four test examples are considered: 

Example 1 [18]. 

y[ = -0.013yi - lOOOyiya, 

y'^ = -2500^22/3, (17) 
y'^ = -0.013yi - lOOOyiya - 2500^2^3, 

t€[0,50], yi(0) = l, 2/2(0) = !, ^3(0) = 0, /lo = 2.9 • 10"^ 

Example 2 [14]. 

y[ = 77.27(^2 - yiy2 + yi- 8.375 • IQ-^y^), 

^2 = (-2/2 - 2/12/2 + y3)/77.27, (18) 
y^ = 0.161(yi-y3), 

t G [0, 300], yi(0) = 4, 1/2(0) = 1.1, ^3(0) = 4, /iq = 2 • 10-^. 



Example 3. 

y[ = -0.04yi + 0.0l2/22/3, 

y'2 = 400yi - 1002/2?/3 - 3OOO2/I, 

2/3 = 30yi, 

t G [0,40], yi(0) = l, y2(0) =y3(0) = 0, /iq = 10"^ 

Example 4. 

y'l = 2/3 - 100yi2/2, 

2/2 = ?/3 + 2?/4-100?/iy2-2-10^yi, 

2/3 = -2/3 + 100yiy2, 

2/4 = -2/4 + lO^yi, 

t G [0,20], yi(0)= 2/2(0) = 1, y3(0) = 2/4(0) = 0, /iq = 2.5 • 10"^ 



The approximation of the Jacobian by a diagonal matrix is used when solving the test problems 
by AS0DE3. For the first test problem the diagonal matrix B with elements 611 = —0.013 — 1000y3, 
622 = — 2500|/3, 633 = — lOOOyi — 25007/2 are used. In the case of diagonal Jacobian approximation 
computational costs of additive methods are dominated by the number of right hand side function 
evaluations. So, computational costs of (2) per integration step are comparable to ones of explicit 
methods. Hence, AS0DE3 is compared with the following numerical codes based on well-known explicit 



RKM4 - 5-stage Merson method of order 4 [15], 
RKF5 - 6-stage Felberg method of order 5 [16], 
RKF7 - 13-stage Felberg method of order 7 [16], 

DPS - 13-stage Dormand and Prince method method of order 8 [17], 
and less well-known Runge-Kutta type method: 
RKN2 - 2-stage method of order 2 [13]. 

The overall computational costs (measured by the number of right hand side function evaluations 
over the integration interval) are given in the table 

Table 1. 

Computational costs of RKM4, RKF5, RKF7, DP8, RKN2, AS0DE3 with stability control. 





Tol 


RKA14 


RKF5 


RKF7 


DPS 


RKN2 


AS0DE3 


1 


10-2 


401 716 


401 005 


982 536 


717 526 


222 441 


243 


10-4 


400 627 


400 656 


982 150 


717 287 


222 481 


5 253 


2 


10-2 


13 391 594 


15 694 434 


38 429 196 


27 998 053 


8 682 849 


4 245 


10-4 


13 384 132 


15 691 105 


38 429 976 


27 993 793 


8 689 861 


89 993 


3 


10-2 


204 889 


237 942 


587 509 


431 591 


133 022 


1 278 


10-4 


206 647 


240 676 


565 396 


430 823 


132 987 


7 908 


4 


10-2 


10 832 


11 874 


29 991 


23 052 


6 585 


174 


10-^ 


10 236 


11 366 


28 819 


23 354 


7 627 


7 938 



8 Conclusions 

In addition to continuum mechanics problems, the constructed additive method can be used for solving 
locally unstable problems. In this case <f{y) corresponds to eigenvalues of the Jacobian matrix with 
positive real parts. In opposite to A-stable methods, explicit Runge Kutta methods are unstable in 
almost the entire right half plane and therefore are more suitable for detecting the local unstable 
solutions. For many locally unstable problems it is also easy to split the right hand side into stiff and 
non-stiff terms from physical considerations. 

So, in this paper, we constructed the third order additive method that is L-stable with respect 
to the implicit part and allows to use an arbitrary approximation of the Jacobian matrix without loss 
of accuracy. Automatic stepsize selection based on local error and stability control are performed and 
the auxiliary formulas for doing this were obtained without significant additional computational costs. 

The aim of numerical computations was to test the reliability and efficiency of the implemented 
integration algorithm with error and stability control as well as with diagonal .Jacobian approximation. 
They didn't aim at solving practical problems of continuum mechanics and locally unstable problems. 
Numerical experiments show reliability and efhciency of the presented method. It follows from them 
that the method has good stability properties for solving mildly stiff problems and that the test 
problems turned out to be rather stiff for the explicit Runge-Kutta methods considered above. It is 
worth remarking that computational costs per step are comparable for both the additive method (with 
diagonal Jacobian approximation) and explicit ones. So, the implemented integration algorithm makes 
it possible to expend the range of applicability of explicit Runge-Kutta methods towards more stiff 
problems. 
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